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Abstract 

We present Monte Carlo simulations of semidilute solutions of long self- 
attracting chain polymers near their Ising type critical point. The polymers 
are modeled as monodisperse self-avoiding walks on the simple cubic lattice 
with attraction between non-bonded nearest neighbors. Chain lengths are up to 
= 2048, system sizes are up to 2"^^ lattice sites and 2.8 x 10^ monomers. These 
simulations used the recently introduced pruned-enriched Rosenbluth method 
which proved extremely efficient, together with a histogram method for estimat- 
ing finite size corrections. Our most clear result is that chains at the critical 
point are Gaussian for large A^, having end-to-end distances R ~ y/N. Also 
the distance Tq — Tc{N) (where Tq = limjv_>oo 7c(A)) scales with the mean 
field exponent, Tq — Tc{N) ~ 1/y/N. The critical density seems to scale with 
a non-trivial exponent similar to that observed in experiments. But we argue 
that this is due to large logarithmic corrections. These corrections are similar to 
the very large corrections to scaling seen in recent analyses of 0-polymers, and 
qualitatively predicted by the field theoretic renormalization group. The only 
serious deviation from this simple global picture concerns the A^-dependence of 
the order parameter amplitudes which disagrees with a minimalistic ansatz of de 
Gennes. But this might be due to problems with finite size scaling. We find that 
the finite size dependence of the density of states P{E, n) (where E is the total 
energy and n is the number of chains) is slightly but significantly different from 
that proposed recently by several authors. 
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1 Introduction 



Consider long flexible polymers in a not too good solvent. At high temperatures, they 
will form extended random coil configurations which can be modeled by self avoiding 
random walks. When T is lowered, there will be a critical temperature where the 
chains start to coagulate and the liquid unmixes. This increases with chain length 
A^. At the theta temperature defined as Tq = limTv^oo Tc{N) also a single (but infinitely 
long) chain will collapse. Qualitatively, this unmixing is described by the mean field 
theory developed by Flory and Huggins The phase diagram is sketched in fig.l. 

For any finite A^ we expect that the internal structure of the dissolved "particles" 
(=chains) becomes irrelevant in the infrared limit, i.e. close to Tc{N). Thus, the 
critical behavior should be fully described by the Ising model, i.e. by the 0{n) sigma 
model with n = 0, |^. This is indeed supported by all available evidence, but it 
tells only half of the story. As in any critical phenomenon there are universal and non- 
universal properties. While the (Ising-)universal aspects (critical exponents and scaling 
functions) must be independent of A^, all non-universal parameters should depend on 
A^ systematically, and can be expected to display another type of scaling in the limit 
A^ —>■ oo. Also, the critical region must become smaller and smaller as A^ is increased, 
and there must be a cross-over to the critical behavior holding at the G point. The 
latter is formally described by a tricritical point in the 0{n) model with n = 0. As any 
tricritical point in 0{n) models, its upper critical dimension is d = 3, whence it should 
display mean field behavior with logarithmic corrections p[. 




Figure 1: Schematic phase diagram for semi-dilute solutions of chain polymers. The up- 
permost curve gives the monomer concentration inside an infinitely large collapsed globule 
under zero outside pressure. The lower curves are coexistence curves for fixed chain length 
N. These curves are strictly monotonically ordered, with the coexistence curve for N below 
that for A^' if iV < iV'. 
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The cross-over should be described by scahng laws which are accessible to experi- 
ment, but which cannot yet be obtained fully from the field theoretic renormalization 
group. Thus even the basic critical exponents are not known, and are the subject of 
controversial speculations. The main theoretical difficulty is that the Ising transition 
and the G-collapse have different upper critical dimensions, making an e-expansion 
non-trivial. 

Let us denote the monomer density by 0. In Flory-Huggins theory one assumes 
that the entropy per unit volume is given by [0 

^=-^ln^-(l-0)ln(l-0), (1) 

and the energy is a quadratic function of independent of A^. Keeping only the mixing 
entropy, expanding ln(l — (p) in powers of (f), and dropping the irrelevant linear term in 
(p, the free energy per unit volume is thus given by 

f3F = ^ln<p+^v<j)^ + ^w4>^ + ... . (2) 

The B point corresponds to a vanishing of v. We assume that f is a linear function 
of T, while > is constant. All (mean field) scaling laws can be obtained from this 
ansatz. In particular, the critical density and the distance from the G point both turn 
out to scale as N~^^'^, 

<pc ~ 1/ViV (3) 

and 

Te-UN)r^l/VN. (4) 
For N = oo, the density at T < Te scales linearly as 

0OO ~ - T. (5) 

Unfortunately, the theory - being mean field - predicts also the classical value 1/2 
for the order parameter exponent f3. For finite A^, the density difference along the 
coexistence curve (the 'order parameter') is predicted to satisfy a scaling law 



Ar"V2/((T,-T)v^) (6) 



with f{x) ~ a/x for x ^ 0, and f{x) ~ x for x — > oo. A minimal modification which 
gives the correct value of P (~ 0.325) was proposed by de Gennes 0]. He suggested to 
keep eqs.(|) to (^, and to adopt simply a different behavior for /(x): 

/(x) ~ x^ , X ^ . (7) 

This ansatz seems theoretically consistent although it cannot be derived from an 
underlying microscopic theory. But unfortunately eq.(^) is in serious conflict with 
experiment. Most experiments 0, H, |, 0, H and subsequent analyses H, |10|, |ll], ^ 
agree that 



0c ~ N~-' (8) 
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with 

= 0.38 ± 0.01 . (9) 
This discrepancy has given rise to a number of speculations and more or less well 



founded conjectures [QQ, 11, 13, 04, Ma, Ra, 17|. One conjecture is that the chains might 



be partly collapsed at the critical point. Since Tc{N) Tq for N ^ oo, the natural 
first assumption is that the end-to-end distance or the radius of gyration should scale 
as R ^ VN, i.e. the chains should be free. One might question this since actually 
Tc{N) < Tq for any finite A^. Accordingly, it was suggested in |]T6[ that chains are 
partly collapsed, 

Rn ~ N""" (10) 

with an exponent xq < 1/2. Actually, this is most unlikely. On the one hand, the 
effective collapse temperature for finite A^ - defined as that T where Rn is the same 
value as for ideal chains - is below Tq [^. On the other hand, chains must be in 
contact at the critical point (otherwise they would not interact) and can penetrate 
each other. Thus there is no force which should compress them beyond the ideal shape 
which maximizes entropy. This argument is very similar to that which explains why 
chains in dense melts are basically ideal. 

In spite of these doubts, it seems wise to leave the exponents eqs.(^ and (|T^) open 
at the present point, and to assume similar scaling laws also for the other observables: 

0(2) _ 0(1) _ - rfN-""' (11) 

and 

TQ-T,{N)r^N-^^ . (12) 

Notice that eqs.(|D and (|^) would give xi = (1 — [3)12 ^ 0.34. Experimentally, the 



exponents are |^ |TB| 

xi ^ 0.23 - 0.34, X3 f« 0.47-0.5. (13) 

In addition, we can generalize eq.(|) to 

0OO ~(Te-T)^ (14) 

with some unknown exponent y, although it seems that all authors have assumed that 
y = 1, with a single exception to be discussed below. 

Since the G-point is a tricritical point, one should expect logarithmic corrections 
to these scaling laws |T^, To leading order, they are obtained by replacing v and 



w in eq.(|^) by their renormalized values for large A^, 

V const (T-Te) [In A^] w const/ \nN. (15) 



Among others this gives [IS 



[In ATI 

^ ^ (16) 



A^ 
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and 

Te - UN) ~ iV-V2[iogiV]-3/ii , (17) 

From this follows 0, ~ (T© - T,)[log(re - T)Y/^\ 

Similar corrections are predicted for the end-to-end distance of single chains at the 
0-point, for the specific heat, and for other observables |Q (but not for pure Ising 
properties such as (f)^"^^ — (f)^^^). These corrections for single chains are not (yet) fully 
seen in simulations. More precisely, simulations [^, ^ show for most observables 
corrections which are much larger than those predicted by theory, with only weak 
hints that the theoretical predictions are correct for extremely large |]22[| . 

One such case where corrections to mean field behavior are very large is the de- 
pendence of 0OO on Tq — T. Based on simulations of very long chains, it was shown in 
|21[] that a best fit is given by y = 0.7. But using even longer chains, with up to 



10^, it was concluded in ||22| that this is an effect of very large corrections to scaling. 



and that y = 1 is indeed the most likely value. Indeed, as shown in |T^, one expects 
logarithmic corrections with the same power as for 0c, 



(re-T)[log(Te-r)] 
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Let us forget for the moment any logarithmic corrections, and assume that the 
scaling laws (H) - (^4|) hold. We should point out that the exponents Xi and y are not 
independent. A simple argument gives 

xsy <X2 . (19) 

To derive this, one just needs that (pc{N) < 0oo at fixed temperature T = Tc{N). But 
this follows from the fact that phase coexistence regions grow with A^: if the point 
(T, 0) is in the phase coexistence region for some A^, it is also in this region for all 
N' > N (we have not found a mathematically rigorous proof for this, but it seems 
heuristically obvious, and it was assumed tacitly by all previous authors). 

Inserting y = I and the numerical values (P) and (H), we see that eq. (|T9|) is violated. 
It seems that this simple observation has been overlooked in all recent literature, and 
it makes several claims obsolete. One either has to admit that y < 1, or that X2 is 
much closer to its mean field value than suggested by experiments, or (which seems 
the least likely) that x^ deviates strongly from its mean field value. 

In view of this unclear situation we decided to perform large scale simulations using 
a novel Monte Carlo algorithm, the Pruned- Enriched Rosenbluth Method (PERM) ||22|| . 
We refer to this reference and to sec. 2 for a description of this algorithm. It allowed us 
to simulate very large systems, with chains of length up to 2000 and beyond. Indeed, 
we could have gone even further as concerns chain length, but we had problems going 
to systems containing more than ^ 400 chains, even if these chains are short. This 
is due to the fact that PERM performs excellently at the 0-point, but becomes less 
efficient at temperatures much below T©. The latter is needed for simulations of short 
chains at the critical point. 
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Thus all our simulations were done on lattices with finite volume \^ in a regime 
where finite size corrections are important. To perform a detailed finite size analysis, 
we used the histogram method as proposed in and applied to polymers in 

|[T8[| . In this method, one constructs the microcanonical distributions^ ("histogram") 
P{E, n) where E is the energy and n is the number of chains. Ising universality requires 
that P{E, n) has certain scaling properties which can be used to extrapolate to the 
thermodynamic limit. In the Ising model proper (which is just the = 1 limit of the 
lattice polymer model studied below), 2nN /V — 1 is replaced by the magnetization M, 
and P{E, M) is symmetric under M — ^ —M. This symmetry is broken for > 1, 
but restored at the critical point (i.e., for E E^ and n «i in the thermodynamic 



limit. It was claimed in [23, 24, 18 1 that the dominant symmetry breaking term in the 



vicinity of (i?c, "^c) can be removed by an affine transformation 

S = E-rn , U = n- sE (20) 



with A^-dependent parameters r and s ("field mixing" |£5l). We found that this is 
not true in our case. Although parameters r, s can be found such that the marginal 
distributions P{M) and P{£) agree with the Ising universal curves within error bars 
(implying also that P{M) is symmetric around A^;), the 2-dimensional distribution does 
not become more symmetric by this transformation. Nevertheless, the possibility to 
compare with the precisely known critical magnetization distribution of the Ising helps 
enormously in fixing the critical point, and extracting critical parameters. We thus 
basically verified that the histogram method gives very reliable (and large!) finite size 
corrections, allowing us to obtain precise critical parameters from fairly small system 
sizes. 

The algorithm and computational details are discussed in the next section. In 
sec. 3 we describe the histogram method, and our main results are presented in sec. 4. 
One interesting aspect of PERM is that it allows to compute free energies with high 
precision. We use this to compare directly with the mean field ansatz eq.(0). We 
conclude with a discussion of our results in sec. 5. 



2 Algorithm and Computational Aspects 



The Pruned-Enriched Rosenbluth Method (PERM) ||2^ is a chain growth algorithm 
based on the Rosenbluth-Rosenbluth (RR) method ||26[ and on "enrichment" or 



copying of successful partial chains. As is well known ||28[, the main drawback of the 
RR method is that it leads to weighted samples with very uneven weights. This is 
counterbalanced in PERM by enrichment (when a large weight chain is copied k times, 
all /c + 1 copies receive l/(/c + 1) times the original weight) and by pruning: low weight 
chains are either pruned (deleted) or, with the same probability, doubled in weight. 



Similar algorithms have been used in [29, 30|. There, however, they were implemented 



"'^Strictly spoken, P{E,n) is the microcanonical partition sum multiplied by e '^'■^ and nor- 
malized so that J2e n P{^^ n) — I. 
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in a 'breadth first' way, by keeping a large population of chains simultaneously in 
memory. Such a strategy would have been not feasible for the very large systems (up 
to several hundred thousand monomers) treated in this paper. Instead, we used a 
recursive 'depth first' strategy. The basic algorithm is described in detail in p2[, some 



further tricks to make it faster and more robust are discussed in [^, 0. We refer to 
these papers for details. We just mention that we hand-coded the recursion without 
using recursive function calls, in contrast to the algorithm shown in the appendix of ||22|| , 



since we otherwise would have had problems with storage. All simulations reported in 
this paper were done on workstations and used less than 50 MB main memory. 

As designed originally, PERM applies to classical (distinguishable) particles which 
are tied together to form a chain polymer. To simulate systems with several polymer 
chains which are indistinguishable and not connected, we need some modifications. The 
first is that each time when we start a new chain, we have to chose the location of its 
first monomer anywhere in the not yet occupied volume. For the k-th chain in volume 
V, this implies a Rosenbluth factor V — {k — 1)N. The total Rosenbluth weight for a 
system consisting of n chains is then nfc<n(^ ^ ^ 1)^) times the product of proper 
Rosenbluth factors for the second, third, ... A^-th monomers in the chains. Secondly, 
we have to multiply the final weight by 1/n! since there are exactly n\ possibilities to 
build up a configuration of n chains sequentially. 

In PERM one uses the product of Rosenbluth and Boltzmann weights of partially 
assembled systems to steer doubling and pruning. In order to have a smooth depen- 
dence of the doubling and pruning thresholds on the number of already assembled 
monomers, we omitted the above factors during the growth, and added them later 
when the system was already built up. 

All simulations were done on the simple cubic lattice with helical boundary con- 
ditions (lattice sites are indexed by a single integer i, and i + V = i). Chains were 
modeled by self-avoiding walks with attractive energy e = —l/ks between each pair 
of neighboring non-bonded monomers. Chain lengths were powers of 2. Volumes were 
also powers of two, V = with m such that systems at the critical point had roughly 
100 chains. Thus V changed from 2^^ for A^ = 8 to 2^1 for A^ = 2048. CPU times (on 
DEC Alpha machines with 400 Mhz) ranged from a few hours for A^ = 8 to roughly 
two weeks for A^ = 2048. 

In these simulations we used one constant chemical potential /i for each chain, 
and another potential /x' for each monomer. Since we measured observables only after 
chains had been finished, this corresponds to an overall fugacity z = e^e^^ per chain. 
We measured: 

• The canonical partition sum Z„ for each n = 1,2,..., nmax- This is needed 
anyhow for the control of doubling and pruning. 

• The average energy En of the total configuration containing n chains, n = 

1,2,..., Tiraax' 

• The average energy e„ of the n-th inserted chain at the moment of its insertion. 
Notice that is approximately equal to J2k=i ^k, but not exactly. 



7 



• The average r.m.s. end-to-end distance RN,n, averaged over all chains. 

• The end-to-end distance rN,n for the last inserted chain. Again, RN,n is approxi- 
mately but not exactly equal to r7v,n- 

• A histogram which contains for each pair [E, n) the microcanonical partition sum 
multiplied by e'^^z^. 

We do not present results for chains shorter than = 8, mainly because this 
requires simulations at rather low temperatures. The algorithm becomes inefficient 
there. For instance, for the Ising model (A^ = 1) we had no problems to simulate 
lattices of size 4^, verifying thereby that the algorithm works in principle also on this 
extreme case. But already for systems with \^ = 8^ we encountered problems. This is 
not too surprising. In spin language, we start simulations with all spins down. We chose 
random sites where spin is still down and flip it up, keeping track of the Boltzmann 
weights by giving each configuration a weight different from 1. If this weight becomes 
too large or too small, we copy or prune, respectively. Finally, if pruning did not 
kill us, we end up at the configuration with all spins up. The total weight of this 
state should be the same as that of the starting configurations with all spins down. 
While this symmetry is exactly respected by the Rosenbluth method without pruning 
or copying, it is not by PERM due to the stochastic nature of pruning. In the long 
run, of course, the symmetry will be approached closer and closer. But this is due to 
very rare events with very large weights. It not only involves huge statistical errors 
but also a systematical bias unless one is very careful and does not estimate statistical 
errors from fluctuations within small samples. These problems have to be kept in mind 
when dealing with A^ > 1. In this case there is no symmetry to check, but there is the 
same danger of underestimating the contribution of rare outlyers. We hope to have 
minimized this danger by using very large samples, typically 10^ to 10^ independent 



"tours" in the terminology of [22 . 



3 Histogram Method and Finite Size Scaling 



Let us consider for the moment an Ising system in d dimensions with lattice size L. We 
denote by Zl{T, h) the canonical partition sum, and by Zl{E, M) the number of states 
with energy E and magnetization density M. The histogram Pl{E,M) is defined as 



PLiE,M) 



-PE+hL'^M 



Zl{T, h) 



-Zl{E,M) . 



For vanishing external magnetic field h and at T = Tc it should scale as ^ 
Pl{E, M) ^ L^/'^-^/'^^(L^/^M, {E - Ec)L'^^'') 



(21) 



(22) 



where E^ is the average energy at the critical point, and (3 = 0.327 ± 0.001 and u = 
0.630 ± 0.001 |3^ are the usual critical indices. The function g is universal up to 
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rescalings of its value and arguments by arbitrary factors. Summing over E resp. M 
we obtain the distributions of magnetization resp. energy, 



Pi(M) ^ L^/'^pMiL^^^M) (23) 

and 

PLiE) ^ L-'/'^PeUE - E,)L-'/n . (24) 

The energy scahng function pe{x) has a single maximum. In contrast, the magnetic 
scaling function pm{x) has two (symmetric) maxima for d < A. 

This is the main point which a finite size scaling analysis has to take into account 
correctly. Naively, the critical point is defined as that temperature where Pm{x) changes 
from being single-humped to double-humped. But strictly, this is true only in the ther- 
modynamic limit, and analyses which neglect this when estimating Tc |3^, ^ can 
have very large systematic errors. More precisely, in c? = 3 the ratio Pm{0) / msix^ pm{x) 
is roughly 0.44, with an error probably < 0.02 In the (i?, M)-plane, the distribu- 
tion Pl{E,M) at h = t = has two peaks located on the two halves of a U-shaped 
support. 

Let us now consider a system like a critical gas where M is replaced by the particle 
density p (= (p/N), and the symmetry M — — M is lost. In the following we shall use 



the total particle number n instead of p. According to p3|, one just has to replace 
E and n by linear combinations 

£ = E-rn (25) 

and 

M = n-sE (26) 

with suitable constants r and s. One arrives at scaling laws for £ and M which formally 
coincide exactly with eqs. (p2| - plD except for the fact that also Mc = (A/") is different 



from zero and that M — Mc-, being an extensive quantity, replaces not M but L'^M: 

Pl{£M) ~ L^'^-^'^-^giiU -K)L^/''~\ {£ - £c)L"^''') . (27) 

This conjecture was tested in ^ |1^] mainly by verifying that the projected 1- 
dimensional distributions Pl(A/') and Pe{£) agreed numerically with the Ising scaling 
functions. For polymer solutions we expect r and s to scale with in a universal way. 



Notice that the justifications for eqs.(^) and (p6[) are rather different. Replacing 
the energy density E' by a linear combination E — rn just corresponds to changing the 
internal energy of the particles. Since that is arbitrary anyhow, we see that eq . (pSD 
is very natural. It just corresponds to a shift in the internal energies such that both 
phases have the same energy density per unit volume. No such interpretation can be 
given for eq.(^). 

Let us assume that we have made the simulations at a temperature near Tc. We 
first obtain rough estimates of the critical fugacity Zc and of r by demanding that both 
peaks in Piin) have the same height, and that the two peaks in PL{£,n) occur at 
the same value of £. A typical result, for A^ = 128, is shown in fig.2a. By summing 
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over all E, we obtain the chain number distribution shown in fig.3a. Notice the very 
asymmetric shape of both distributions. Notice also that the ratio between the height 
of the central minimum of Piin) and its maximal value is not exactly 0.44, indicating 
that we have not simulated exactly at T^. 

In a first refinement we now determine a crude estimate of Tc by re-weighting 
PL{S,n) such that the minimum in the projected distribution Piin) has a height 
equal to 0.44 times the average value of the maxima. We then fit s such that PiiAf) 
becomes symmetric under JV — Afc — (A/" — A^), where f/c is the minimum position, 
and readjust at the same time the fugacity Zc such that both peaks have again the 
same height. The fact that we always do find a value of s such that Pl{M) becomes 
approximately symmetric and equal to the Ising scaling function is highly non-trivial. 

The last readjustment will in general have shifted the peaks in the A/')-plane 
such that they occur no longer at the same value of S. We thus repeat the fit of r and 
estimate more precise values of Tc, Zc, and s by repeating the second step. In principle, 
we could iterate this procedure until we reach convergence, but already after the second 
refinement the two peaks of Pl{£,M) will occur at the same £. We then check that 



Pl{M) and Pl{£) agree with the scaling functions given in ||2^ within the statistical 
errors. The final results for the raw data shown in figs. 2a and 3a are shown in figs. 2b 
and 3b. For other chain lengths, results are very similar. 

We thus have succeeded in finding field mixings such that the 1-d projections of the 



2-d histogram agree with the Ising case. Exactly this was done also in p^, |J, |T8|. But 
we see from fig. 2b that the symmetry of the projected density Pl{M) is misleading. 
The 2-d histogram has not become more symmetric by replacing n by M . We have not 
tried to apply formal measures of symmetry to fig. 2, but it seems that fig.2b is rather 
less symmetric than fig. 2a. 

We see thus that linear field mixing is questionable. We tried a number of alterna- 
tives. By far the most successful numerically is the following. We first subtracted from 
E a term linear in n (corresponding just to a redefinition of internal energy), and then 
applied a conformal transformation with a scale factor which depends only on n. We 
did not try to optimize systematically since this is heuristic anyway, but good results 
were obtained with power law factors, 

J\f={n/{n)Yn, £ = {n/ {n)Y{E - r n) . (28) 

The same data shown already in figs. 2a and 2b are plotted against these variables 
in fig.2c, with a = —0.3. We see indeed a much more symmetric distribution. One- 
dimensional projections from this distribution are not very different from those obtained 
with linear mixing, and are not shown. 

These different mixing ansatzes showed us that the estimates of T^, and of the 
critical fugacity are very robust. Heuristically, this can be explained by the fact that 
the support of PL{E,n) is a very narrow band, i.e. the variance of E for fixed n is 
rather small. We thus conclude somewhat paradoxically that we cannot verify a basic 



assumption of p3|, p3, IT8[, but we agree that the histogram method is extremely helpful 



for extracting critical parameters. 
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Figure 2: (a) Contour plot of 
PL{£,n) for chains of length 
A'^ = 256 in a lattice with size 
L = 64. Notice the asym- 
metric shape. For the Ising 
model, the contour plots would 
be symmetric under reflection 
on the vertical axis. The nor- 
malization is arbitrary. (b) 
Same data as in panel a, but 
after transforming n ^ AA by 



means of eq.(26) and readjust- 
ing the other parameters. The 
parameter s is 0.003. (c) Same 
data again, but using the non- 
linear transformation eq.(|2l 
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To verify the correctness of the approach (and of the simulations!), we checked the 
scaling with L for one chain length, = 128. According to eq.(^), the values of 
L~f^/'^ Pi[M) should coincide when plotted against L^I^M . This scaling plot is shown 
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Figure 3: (a) Distribution Piin) = z'^Zn (up to an arbitrary normalization factor) against n, 
obtained by projecting the density of fig.2 onto the horizontal axis, (b) Same as panel a, but 
after transforming n M. To demonstrate the symmetry of the distribution, both curves 
Pii-N") and Pl{2Mc — M) are superimposed. As in the previous figures, the normalization is 
arbitrary. 

in fig.4a, while the analogous scaling plot for the energy distribution is shown in fig.4b. 
We see perfect agreement. For other chain lengths we usually simulated at two lattice 
sizes in order to check for obvious inconsistencies, but we did not check systematically 
the finite size scaling. 



4 Results 



As a first result we show the swelling factor R%/N for three different temperatures as 
a function of the number n of chains (fig. 5). These results were obtained at fixed L 
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and N. The critical point corresponds to the central curve and to ric ~ 79. From this 
we can make several interesting observations: 

1) All chains are swollen, i.e. > VN in all cases. 

2) The swelling is weakest for n = 1, i.e. for isolated chains, and increases as the 
chain density is increased. For large n it seems to saturate at a value < 2. This agrees 
with the fact that all three temperatures are below the G point, and chains at Tq are 
slightly swollen, with R%/N — > 1.7 — 2 for N ^ oo |2l], For large n the chains 
feel mainly the slight imbalance between self avoidance and attraction for short chains 
which makes them swollen like very long B chains, but for small n the attraction from 
the other chains is missing and the chains are slightly less swollen. 

3) There is rather weak dependence on T, and this dependence weakens further as n 
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Figure 4: (a) Scaling plot of the distribution of transformed chain numbers M for N = 128 
and T = 3.168. Within statistical errors, the latter is our estimate for T^- The critical 
monomer density was assume to be (^c = 0.1544. (b) The analogous distribution of the 
transformed energy £. 
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increases. 

As a result of the last point, the systematic uncertainty in R%/N due to the uncer- 
tainty of Tc is negligible. Not negligible, in contrast, is the dependence on the estimate 
of ric- Indeed, due to the very small statistical errors in fig.5, the uncertainty of is 
the largest source of error for .Rtv- 

Results for R%/N at the critical point obtained from fig.5 and from similar figures 
for other values of are shown in fig.6. We see that the swelling increases slightly 



with A^, as also found in []T8|. But in contrast to these authors we see that the swelling 
saturates, ruling out a growth of -Rat with a power > 1/2. Of course, such a growth 
would seem very strange and was rejected on this ground in [l^. But we can even 
extrapolate to A^ = oo (dashed curve in fig.6), and verify that we obtain in this limit 
the same swelling as for isolated G polymers. 

To estimate the exponent xi directly, one would have to simulate much larger 
systems than we could afford. But one can obtain this exponent indirectly from finite- 
size scaling exactly at Tc. Let us denote by AS and AJV the rms widths of Pl{S) 
and pl{J^), respectively. Keeping the A^-dependent non-universal factors omitted in 
eg . (|27D , one easily shows that 



(2) _ ^(1) 



Lj 



(29) 



We tried to use this for estimating Xi. But the results were rather disappointing, 
presumably because of the uncertainty in the field mixing. In the linear mixing ansatz, 
eqs. ( p5tp6D were actually written in p3|, |1^ asamap {S,Af) = {l—rs)~^{E—rn,n—sE). 
We have omitted the factor (1 — rs)^^ up to now since it is irrelevant for everything 
we did up to now, but it does become relevant in eg. (|29|) . But we also had argued that 
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Figure 5: Swelling factors R%/N for fixed = 128 at three temperatures, T = 3.03, 3.1655, 
and 3.30 (bottom to top), plotted against the number of chains n. The central value is close 
to the estimated Tc. All data were obtained with = 2^^ sites. 
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Figure 7: Log- log plot of {Tq —Tc{N)) /Tc{N) against (crosses). The dashed hue has slope 
-0.51. 

eg . (pSD is very natural as it is written, without the factor (1 — rs)~^. Thus we propose 
that the correct linear mixing for our purpose is 



S = E — r n. 



A/" = (1 - rs) ^{n - s E) =n - Y 



-£. 



rs 



(30) 



Using this, we obtain xi ~ 0.18, albeit with large error bars which are hard to estimate 
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since they are a mixture of statistical and systematic errors. 

Alternatively, we also tried the nonlinear transformation eg. (pS]) , and we tried to 
use the untransformed width An instead of AJV. Although the detailed numbers 
entering into eq.(|29|) are then quite different from those obtained with linear mixing, 
the resulting values of xi were again ^ 0.15 to 0.2, confirming thus the estimate 
Xi ~ 0.18. In view of its uncertainty, this estimate is presumably compatible with the 
lower end xi ~ 0.23 of phenomenological analyses, but it seems very hard to reconcile 
it with the prediction xi = (1 - /3)/2 = 0.34 of eqs.(§J3). 

Estimates of Tc{N) are shown in fig. 7. More precisely, this is a log- log plot showing 
{Tq — Tc{N))/Tc{N) against N. We see a very clean scaling law with exponent x^ = 
0.51 ± 0.01, suggesting that the mean field value 1/2 is indeed correct. We should 
however realize that we are not yet very far in the scaling regime, in spite of the large 
values of N. Thus, using slightly different scaling variables as, e.g., (Te — Tc{N))/Tq 
instead of the variable used in fig. 7, might lead to slightly different exponents. Also, 
this plot is very sensitive to the exact value of Tq. In order to reduce this uncertainty 
further, we have made additional simulations of single very long polymer chains using 
the same routine as in [^, but with even larger (up to 1.6 x 10®) and even larger 
system sizes (up to 512^ sites). Our best estimate is now Tq = 3.717 ± 0.002. This 
uncertainty contributes less than 0.01 to the uncertainty of x^. Finally, according to 
eq.(p!7[) we should see weak logarithmic corrections. We believe that we do not see 
them in fig. 7 because of the just mentioned uncertainties. In summary, we can say 
that X3 = 1/2 is the most likely value. 

Finally, in fig. 8 we show our estimates of the critical density. For our largest values 
of A^ (which agree roughly with the longest chains used in experiments) we see a 
power law with exponent X2 ~ 0.38. This agrees perfectly with experiment and with 



phenomenological analyses |3l-[Tl|. If we accept this as the true critical exponent, and 
accept at the same time the mean field exponent y = 1, we run into the problem posed 
by inequality (plQ]). But we see also from fig. 8 that there are very large corrections to 
scaling. If we define effective A^-dependent exponents by fitting locally, they increase 
slightly but systematically with A^. 

This suggests strongly that the deviation from the mean field exponent X2 = 1/2 
is entirely due to finite- A^ corrections which vanish for A^ — 00. This is completely 
consistent with the very large non-asymptotic corrections seen for single chains in the 
limit A^ — >■ cx) in ETI, p2| . To stress the similarity between the critical monomer density 



4>c{N) as a function of Tc{N), and the infinite chain density 0oo at the same value of the 
temperature, we plot them both in fig.9. We see nearly parallel curves which suggests 
that indeed both densities scale with the same power of Te — T. For the same values 
of T, Flory-Huggins theory predicts (poo/ 4>c = 3. For the longest chains our data give 
3.2±0.2 for this ratio, with a slight tendency to increase with A^. We see the same small 
curvature in both curves, suggesting that a pure power fit might not be appropriate. 
Such a fit would give an exponent 0.75 to 0.85, depending on the interval used for the 
fit. The dashed line indicates in contrast the prediction of eg . (|T8D . It does not give 
a perfect fit, but it definitely shows the correct trend. In particular, fitting this curve 
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Figure 8: Log-log plot of (j)c{N) against N. The dashed line has slope 
data for large A^, but there are very substantial deviations at small N. 
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Figure 9: Double logarithmic plots of the infinite- monomer density (poo and of the critical 
densities for finite A^, as functions of Te — T resp. Tq — Tc{N). The values for are mostly 
from |2^, except for the points very close to Tq. The slight scatter of these points reflects the 
dependence on system size. The dashed line corresponds to eq.(|l8|), with r = {Tq — T)/T, 
and with logr arbitrarily replaced by log 1.6r. 



by a pure power law would give an exponent ~ 0.8 to 0.9 (again depending on the 
fit interval), while the correct power is 1. It seems thus very likely that all deviations 
from mean field behavior seen in figs. 8 and 9 are due to logarithmic corrections. 

It seems thus likely that Flory-Huggins theory provides a much better description 
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of unmixing of long polymer chains than previously thought. To check this directly, 
we tested the ansatz eq.(|]) directly. According to it, the free energy for fixed volume 
and temperature should consist of a term which depends only on the monomer density 
0, plus a known entropy contribution which also depends on the chain length N. To 
test this, we made simulations at the same T (= 3.5631) and in the same volume (128'^ 
sites) for chain lengths N = 1024, 2048, and 4096. Within statistical errors, this T is 
the critical temperature for = 2048, thus the chains with = 1024 are deeply in 
the single-phase domain, while those for A^ = 4096 are deeply in the two-phase region. 
This is illustrated in fig. 10a. There, constants are added arbitrarily, and fugacities are 
adjusted arbitrarily such that the curves are flattest in a qualitative sense. 

In fig. 10b we show the same data, but after removing the supposed entropic con- 
tribution —N~^V(f)\og(f). In this panel, additive constants and fugacities are adjusted 
such that the curves coincide for large 0. It is there where eq.(0) should be most reli- 
able: for large monomer densities, where chains penetrate substantially, there should 
be hardly any difference between one chain of length A^ and two chains of length A^/2, 
up to the entropic difference which is taken out in fig. 10b. The non-trivial hypothesis 
underlying eq.(|^) is that the same is true also for small densities. We see from fig. 10b 
that it is not perfectly true, but the A^-dependence in fig. 10b is much weaker than that 
in fig. 10a. A detailed fit shows in addition that none of the curves in fig. 10b can be 
fitted perfectly by a cubic polynomial, showing that the internal energy contains also 
terms ~ 0^ and, since this term has the wrong sign, higher powers. 

Thus the Flory-Huggins ansatz is not exact, but it seems to be a good first ap- 
proximation. We should expect deviations from Flory-Huggins due to logarithmic cor- 
rections, as discussed in sec.l. As pointed out there, the leading corrections preserve 
eq.(2), whence also (poo = ^(pc should still hold. But the coefficients of the quadratic 
and cubic terms, v and w, should decrease slowly with A^. This is indeed found when 
making power law fits to the curves in fig. 10b (both decrease roughly by 20% when go- 
ing from A^ = 1024 to 4096), but we cannot make a more detailed comparison because 
of the presence of higher than cubic terms. 



5 Discussion 

We have applied a novel Monte Carlo scheme to simulate very large systems of chain 
polymers in semidilute solutions. Our system sizes are comparable to those of previous 
analyses as far as chain numbers are concerned. But our chain lengths are much longer, 
extending to > 2000. The latter would have been unfeasible for other algorithms we 
are aware of, and is possible only since our algorithm can make use of the fact that 
long chains close to the critical point are nearly free. 

Our most solid result is that chains at the critical unmixing point are not shrunk. 
They are slightly expanded, but the expansion factor tends to a constant for chain 
length N ^ oo. Thus asymptotically, for N ^ oo, chains are Gaussian in contrast 
to recent speculations, but in agreement with the most recent simulations A 



somewhat less strong result which, however, seems also very clear cut, is that the A^- 
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dependence of the critical temperatures is as predicted by Flory-Huggins (mean field) 
theory. Again this is in contrast to recent speculations. 

Strong deviations from mean field behavior were seen in the critical density 0c- 
Here, a scaling fit would produce the same anomalous critical exponent as seen also 
in experiment. But we show that 0c can obey this seen anomalous scaling only if the 
same anomalous scaling governs also the density inside a very large globule, i.e. a single 
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Figure 10: (a) Total negative free energy —(3FV of systems with three different chain lengths 
(upper curve: N = 1024; middle curve: N = 2048; lower curve: N = 4096) plotted against 
the monomer concentration (p. In all three cases, lattice size and temperature were the same: 
V = 2^^ sites, and T = 3.5631. The latter is close to Tc for N = 2048, while chains with 
N = 1024 (4096) are in the single phase (coexistence) domains, (b) ((p/Nlog(p — PF)V 
for the same systems as in panel a. According to the Flory-Huggins ansatz, this should be 
independent of N. Fugacities and additive constants are fixed such that the curves coincide 
for large (f>. 
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collapsed polymer chain close to the point. For the latter, a superficial analysis also 
suggests anomalous scaling. But theoretical prejudices and more careful simulations 
suggest these might be fake and due to large logarithmic corrections to scaling. These 
corrections are in qualitative agreement with field theoretic predictions, and should 
vanish at extremely large chain lengths and extremely close to the point. 

Our simulations do not suggest that deviations from mean field behavior — which 
must be present because of the anomalous Ising exponents — are "minimal" in the 
sense of de Gennes in contrast to the prediction xi = (1 — l3)/2 = 0.34 (where 
Xi the exponent for the A^-dependence of the order parameter) we find x ~ 0.18, even 
lower than most phenomenological estimates. But we should say that this estimate is 
by far the most shaky of all our results. 

It is not clear whether our predictions can be tested by new experiments or by 
re-analyzing old ones. Our chain lengths are comparable to those in real experiments 

1000 Kuhnian lengths), and we predict that mean field behavior should be seen 
in (pc only for much longer chains. This seems not feasible at present. On the other 
hand, fits involving the logarithmic corrections might show their presence already for 
shorter chains. Measurements of Tc{N) and of chain dimensions could be improved, 
and should be in agreement with mean field behavior since logarithmic corrections seem 
to be small for them. 

The biggest problem in simulations are the finite sizes of the system. In contrast to 
chain lengths, chain numbers in our simulations are orders of magnitude smaller than 
those in real experiments. Nevertheless we believe that finite size effects do not seriously 
affect our conclusions. This is due to the use of sophisticated histogram methods 
which allow a detailed comparison with the finite size behavior of the Ising model. 
Essentially they make scaling ansatzes for the microcanonical partition sum (a similar 
method, but without the correct finite size dependence of the microcanonical partition 
sum, was proposed in [^). We verified this dependence partially, which means in 
particular that we also measured the Ising exponents jS and u within our simulations. 
But we found that the simple linear field mixing proposed in PSI 113] does not work 



particularly well. We showed that a nonlinear mixing ansatz works much better, but 
we have no good theoretical reason for this ansatz. This is an interesting problem 
which deserves further investigations. But it is not very important as far as the precise 
location of the critical point and the extraction of critical parameters are concerned 
(except for the exponent Xi), and it cannot affect the above conclusions. 
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